31  Teleconnections I

Learning Objectives

After completing this tutorial you should be able to

Download the directory for this project here, make sure the directory is unzipped and move it to your bi328 directory. You can open the Rproj for this module either by double clicking on it which will launch Rstudio or by opening Rstudio and then using File > Open Project or by clicking on the Rproject icon in the top right of your program window and selecting Open Project.

There should be a file named 31_elnino_I.qmd in that project directory. Use that file to work through this tutorial - you will hand in your rendered (“knitted”) quarto file as your homework assignment. So, first thing in the YAML header, change the author to your name. You will use this quarto document to record your answers. Remember to use comments to annotate your code; at minimum you should have one comment per code set1 you may of course add as many comments as you need to be able to recall what you did. Similarly, take notes in the document as we discuss discussion/reflection questions but make sure that you go back and clean them up for “public consumption”.

  • 1 You should do this whether you are adding code yourself or using code from our manual, even if it isn’t commented in the manual… especially when the code is already included for you, add comments to describe how the function works/what it does as we introduce it during the participatory coding session so you can refer back to it.

  • Before you get started, let’s make sure to read in all the packages we will need.

    # load libraries ----
    
    # reporting
    library(knitr)
    
    # visualization
    library(ggplot2)
    library(ggthemes)
    library(patchwork)
    library(ggrepel)
    
    # data wrangling
    library(dplyr)
    library(tidyr)
    library(readr)
    library(skimr)
    library(janitor)
    library(magrittr)
    library(lubridate)
    
    # lake models
    library(GLMr)
    library(glmtools)
    
    # mapping
    library(sf)
    library(ggspatial)
    library(rnaturalearth)
    library(rnaturalearthhires)
    library(rnaturalearthdata)
    
    # set other options ----
    options(scipen=999)
    
    knitr::opts_chunk$set(
      tidy = FALSE, 
      message = FALSE,
        warning = FALSE)

    31.1 Teleconnections

    We have already seen how local and regional drivers can impact the thermal structure of lakes. Local impacts could include something like land cover change which alters the stream inflow/outflow of a lake impacting temperature. At a regional scale climate change can alter weather patterns, including both extreme events which disrupt regular seasonal patterns or also gradual change which raises lake temperatures overall.

    Next, we are going to explore how ecosystems can be influenced on a global scale by teleconnections, i.e. meteorological, societal, and/or ecological phenomenon that link remote regions via cause and effect relationships. Simulations are especially powerful in this context, because it allows researchers to manipulate drivers independently from each other to disentangle how they impact ecosystems. By running different teleconnection scenarios, scientists are then able to better predict how ecosystems will respond to remote drivers.

    In our module, we will specifically look at lake thermal structure and ice cover and determine how it is impacted by climate teleconnections to decadal events such as the El Niño/Southern Oscillation (ENSO). This is an example of a global driver. ENSO period with warmer ocean surface temperatures in the Pacific Ocean affect regular trade wind patterns. This altered atmospheric circulating impacts air temperature and precipitation patterns globally. How ENSO ends up impacting a particular lake will be further mediated by local and regional dynamics.

    Our central question is:

    How do global ENSO teleconnections interact with local & regional drivers to affect lake temperatures, thermal structure, and ice cover?

    We can address this question using the General Lake Model we used to understand the effect of local & regional drivers on the thermal structure of lakes. To explore those interactions, we varies the local/regional drivers (weather/climate) and ran different scenarios (described in our met file) but kept the lake we were using constant, i.e. the morphology of the lake and other specific characteristics described in the nml file did not change.

    ENSO has large-scale effects, however, the way an El Nino year plays out will vary by region, additional local effect will further modify the impacts on the thermal structure of each lake.

    Consider this

    Regional differences in effect of El Nino

    During El Nino years there is a shift in the climate pattern of the Pacific Ocean, resulting in unusually warm surface waters in the the eastern equatorial Pacific Ocean. Use the figure above to describe regional differences in climate patterns in North America are impacted.

    Did it!

    [Your answer here]

    To be able to explore our question we will need nml files describing lakes in different regions across North America. You should a folder called sim_lakes to your project directory. If you take a look in there, you fill find individual folders for a set of lakes, in each folder you have an nml file describing that lake along with a met file containing observed meteorological data. These are all lakes that are part of the GLEON network.

    Let’s take a look at where those lakes are situated.

    lakes <- read_delim("data/lake_characteristics.txt", delim = "\t") %>%
      clean_names()
    
    kable(lakes)
    lake_name state latitude longitude elevation_m max_depth_m surface_area_km2 trophic_state mixing_regime watershed_land_use neon_domain data_provider links
    Barco Lake Florida 29.67594 -82.00823 28 7 0.12 Oligotrophic Warm monomictic Forest 3 NEON Core Aquatic http://www.neonscience.org/field-sites/field-sites-map/BARC
    Crampton Lake Wisconsin 46.21025 -89.47345 591 18 0.26 Oligotrophic Dimictic Forest 5 NEON Core Aquatic http://www.neonscience.org/field-sites/field-sites-map/CRAM
    Falling Creek Reservoir Virginia 37.30333 -79.83722 507 9 0.12 Eutrophic Dimictic Forest 7 GLEON; Carey Lab https://www.westernvawater.org/drinking-water/water-sources-and-treatment
    Lake Mendota Wisconsin 43.10970 -89.42060 259 25 39.39 Eutrophic Dimictic Agricultural/urban 6 LTER- North Temperate Lakes https://lter.limnology.wisc.edu/researchsite/lake-mendota
    Lake Sunapee New Hampshire 43.38020 -72.05270 333 33 16.74 Oligotrophic Dimictic Forest 1 GLEON; Lake Sunapee Protective Association http://www.lakesunapee.org/lake-sunapee-history/
    Prairie Pothole North Dakota 47.13012 -99.25302 565 3 0.11 Mesotrophic Dimictic Grasslands 9 NEON Core Aquatic http://www.neonscience.org/field-sites/field-sites-map/PRPO
    Suggs Lake Florida 29.68770 -82.01785 29 6 0.73 Mesotrophic Warm monomictic Forest 3 NEON Core Aquatic http://www.neonscience.org/field-sites/field-sites-map/SUGG
    Toolik Lake Alaska 68.62956 -149.61051 738 26 1.46 Oligotrophic Dimictic Tundra 18 NEON Relocatable Aquatic; LTER- Arctic http://arc-lter.ecosystems.mbl.edu/toolik-main
    Table 31.1: Characteristics for eight lakes in the GLEON network.

    Let’s make a map displaying their geographic locations. We’ve previously used geom_text() to add text (labels) to maps. But those weren’t exactly pretty maps. We are going to add a new tool to our toolbox for generating maps, geom_label_repel(). This function works especially well if you want to combine using a shape to plot a point and add a label that is placed in a way that it does not overlap with that point or other labels on the plot.

    # dowhload world map as shapefile
    world <- ne_countries(scale = "medium", returnclass = "sf")
    
    # download us states as suhape file
    us_states <- ne_states(country = "United States of America", returnclass = "sf")
    
    # lat/long for map extent
    x_min <- min(lakes$longitude) - 2
    x_max <- max(lakes$longitude) + 2
    y_min <- min(lakes$latitude) - 2
    y_max <- max(lakes$latitude) + 2
    
    # set color for fill
    map_color <- "khaki3"
    
    # create plot
    ggplot() +
      geom_sf(data = world, color = "black", fill = map_color) +  # plot outline of countries
      geom_sf(data = us_states, color = "black", fill = NA) +     # plot outline of states
      coord_sf(xlim = c(x_min, x_max),
               ylim = c(y_min, y_max)) +                          # plot boundaries for map
      geom_point(data = lakes, aes(x = longitude, y = latitude),  # add sites
                size = 2) +
      geom_label_repel(data = lakes, 
                       aes(x = longitude, y = latitude, label = lake_name),
                       size = 3)

    Lake locations
    Consider this

    Use the map describing changes in regional climate patterns due to El Nino and the map of the distribution of Lakes to hypothesize how you think lakes will be affected during El Nino years compared to normal years. Consider both how lake temperatures and ice cover could be impacted. Be specific - your answer should consider (among other things) which lake you would expect to display the strongest/weakest degree of warming during an El Nino year, which lakes might not be affected. Consider that we have already discovered that patterns of precipitation and wind also impact thermal structure of lakes.

    Did it!

    [Your answer here]

    Consider this

    Briefly describe how you can design an experiment to determine whether El Nino impacts the thermal structure of lakes and identify the effect of regional/local drivers modifying that global interaction. Be specific - consider what information (data sets) you need and what simulations you would run to compare.

    Did it!

    [Your answer here]

    Don’t worry - you won’t have to run simulations for every lake - we will divide and conquer!

    31.2 Set up simulation

    If you take a look at the sim_lakes folder, you should see that it contains individual folders for a series of lakes. In this example, we will use Lake Sunapee - this means that when you run the code for your lake, you will need to be sure to change the lake name accordingly.

    First, let’s create a vector with our lake name2.

  • 2 The advantage for doing it this way instead of directly in the individual code lines is that we only have to set it once, this makes the code a lot more reusable because we don’t have to dig through the code and find all the arguments that set a file path or specify the lake name.

  • LakeName <- "Sunapee"

    Next, we will set the file path for the sim folder. Each lake sim folder contains the input files needed to run GLM: the nml file which describes the lake, the met file describing the drivers of thermal structure, and input files describing inflow and outflow to the lake.

    sim_folder <- file.path("sim_lakes", LakeName)

    Next, we’ll specify the file path for the nml file, read it in, and take a look at it.

    # define file path
    nml_file <- file.path(sim_folder, "glm2.nml")
    
    # read nml file
    nml <- read_nml(nml_file)
    
    # view file
    print(nml)
    &glm_setup
       sim_name = 'MacrosystemsEDDIE_Sunapee'
       max_layers = 250
       min_layer_vol = 0.025
       min_layer_thick = 0.05
       max_layer_thick = 0.75
       Kw = 0.17
       coef_mix_conv = 0.08213
       coef_wind_stir = 1.27
       coef_mix_shear = 0.296
       coef_mix_turb = 0.526
       coef_mix_KH = 0.415
       coef_mix_hyp = 0.453
    /
    &morphometry
       lake_name = 'Sunapee'
       latitude = 43.395
       longitude = -72.053
       bsn_len = 12870
       bsn_wid = 3000
       bsn_vals = 82
       H = 299.43, 299.943, 300.443, 300.943, 301.443, 301.943, 302.443, 302.943, 303.443, 303.943, 304.443, 304.943, 305.443, 305.943, 306.443, 306.943, 307.443, 307.943, 308.443, 308.943, 309.443, 309.943, 310.443, 310.943, 311.443, 311.943, 312.443, 312.943, 313.443, 313.943, 314.443, 314.943, 315.443, 315.943, 316.443, 316.943, 317.443, 317.943, 318.443, 318.943, 319.443, 319.943, 320.443, 320.943, 321.443, 321.943, 322.443, 322.943, 323.443, 323.943, 324.443, 324.943, 325.443, 325.943, 326.443, 326.943, 327.443, 327.943, 328.443, 328.943, 329.443, 329.943, 330.443, 330.943, 331.443, 331.943, 332.343, 332.443, 332.543, 332.643, 332.743, 332.843, 332.943, 333.043, 333.143, 333.243, 333.343, 333.443, 333.543, 333.643, 333.743, 333.943
       A = 1, 16.90309, 38.87712, 64.23176, 158.88908, 7254.80801, 25205.89398, 33207.81875, 41551.18602, 51853.6219, 62734.1436, 79030.41666, 105037.51731, 133941.80829, 168576.24818, 214836.63623, 269137.82616, 327112.05845, 401221.9844, 483528.22071, 575947.57823, 676930.04325, 797467.69789, 925704.7119, 1088626.8751, 1283336.99713, 1505854.39934, 1740820.93038, 1988527.32348, 2234911.89462, 2485285.59597, 2723309.89675, 2978076.71227, 3264022.97524, 3653258.97559, 4050485.06854, 4414109.50064, 4796077.1705, 5158332.45198, 5579849.98152, 6016348.72339, 6468883.43068, 6928267.27171, 7476296.00935, 8028820.97004, 8623779.45825, 9274820.72252, 9852960.63206, 10365149.7392, 10843023.87491, 11253594.96068, 11644117.28649, 12013032.38706, 12356310.56475, 12673823.35603, 12988489.66626, 13285122.06535, 13576406.32546, 13863498.61821, 14148319.13512, 14418854.84716, 14691317.51194, 14951571.07182, 15234317.57908, 15510517.51797, 15780918.00526, 16009460, 16071489, 16135107, 16203453, 16273655, 16349525, 16489737, 16603406, 16760708, 16788809, 16804233, 16826128, 16844170, 16863779, 16885473, 16934251.6
    /
    &time
       timefmt = 2
       start = '2012-11-01 00:00:00'
       stop = '2013-12-31 23:00:00'
       dt = 3600
       timezone = -5
    /
    &output
       out_dir = '.'
       out_fn = 'output'
       nsave = 24
    /
    &init_profiles
       lake_depth = 33.7
       num_depths = 6
       the_depths = 0, 4, 8, 12, 16, 20
       the_temps = 11, 11, 11, 11, 11, 11
       the_sals = 0, 0, 0, 0, 0, 0
    /
    &meteorology
       met_sw = .true.
       lw_type = 'LW_IN'
       rain_sw = .false.
       atm_stab = .false.
       catchrain = .false.
       rad_mode = 2
       albedo_mode = 1
       cloud_mode = 4
       meteo_fl = 'met_hourly.csv'
       subdaily = .true.
       wind_factor = 1
       sw_factor = 1
       lw_factor = 1
       cd = 0.0013
       ce = 0.0013
       ch = 0.0013
       rain_threshold = 0.01
       runoff_coef = 0.3
    /
    &inflow
       num_inflows = 1
       names_of_strms = 'stream'
       strm_hf_angle = 65
       strmbd_slope = 3
       strmbd_drag = 0.016
       inflow_factor = 1
       inflow_fl = 'inflow.csv'
       inflow_varnum = 3
       inflow_vars = 'FLOW','SALT','TEMP'
    /
    &outflow
       num_outlet = 1
       flt_off_sw = .false.
       outl_elvs = 331
       outflow_fl = 'outflow.csv'
       outflow_factor = 1.05
       seepage_rate = 0
    /
    &snowice
       snow_albedo_factor = 1
       snow_rho_max = 500
       snow_rho_min = 100
    /
    &sed_heat
       sed_temp_mean = 12
       sed_temp_amplitude = 11
       sed_temp_peak_doy = 242.5
    /

    One of the arguments in the nml file specifies the file with the information on drivers (i.e. it directs to the met_hourly.csv) file. Let’s take a look at the meteorological input data for the duration of the simulation.

    plot_meteo(nml_file) 

    Figure 31.1: Metereological data set summarizing the behavior of drivres of thermal structure for the baseline simulation.
    Consider this

    Briefly describe the meteorological drivers over the course of the simulation. Be specific (what are ranges for each value? How do they change over the course of the year?)! Consider whether those make sense given the location of your lake (think biomes/climate zones).

    Did it!

    [Your answer here]

    31.3 Run & analyze baseline simulation

    Now that we have everything set up, we can run the model.

    run_glm(sim_folder, verbose = TRUE)
    [1] 0

    To analyze the output, let’s specify the file path for the output file that is now in our lake simulation folder (output.nc). Again, this will be our baseline data set based on observed meteorological data for a specific year (in our example 2013).

    baseline <- file.path(sim_folder, "output.nc")

    Let’s take a look at the variables that were output as part of our GLM simulation run.

    # pull variable names from output
    var_names <- sim_vars(baseline)
    
    # print variables in output
    kable(var_names,
          caption = "Variable names and units for lake simulation.") 
    name unit longname
    evap m/s evaporation
    extc_coef unknown extc_coef
    hice meters Height of Ice
    hsnow meters Height of Snow
    hwice meters Height of WhiteIce
    I_0 10E-6m Shortwave
    NS Number of Layers
    precip m/s precipitation
    rad unknown solar radiation
    rho unknown density
    salt g/kg salinity
    temp celsius temperature
    Tot_V m3 lake volume
    u_mean m/s mean velocity
    u_orb m/s orbital velocity
    V m3 layer volume
    wind m/s wind
    Variable names and units for lake simulation.

    One of the variables we are especially interested in is hice which describes the lake thickness on the surface of the lake.

    Part of our central question is understanding how water temperatures at the lake surface and bottom (below the thermocline) be affected by El Nino teleconnections.

    To estimate the bottom depth we can extract the lake depth in meters from the output.

    # extract lake depth
    LakeDepth <- get_surface_height(baseline)
    
    head(LakeDepth)
        DateTime surface_height
    1 2012-11-02       33.70000
    2 2012-11-03       33.69793
    3 2012-11-04       33.69418
    4 2012-11-05       33.69017
    5 2012-11-06       33.68645
    6 2012-11-07       33.68364
    Give it a whirl

    Create a table with summary descriptive stats for the lake depth (min, max, mean, std) and describe the level of observed fluctuation over the course of the simulation.

    Did it!

    [Your answer here]

    Here’s what your table could look like.

    mean_depth std_depth min_depth max_depth
    33.7 0.2 33.4 33.9
    Table 31.2: Descriptive statistics of lake depth for Lake Sunapee, NH

    We are going to want to extract the lake temperature data from our baseline simulation into a data.frame so we can compare it to our teleconnection scenarios down the line. We will extract it at the surface and for the bottom by setting the argument z_out to 0 (surface) and the minimum observed depth over the simulation (min(LakeDepth$surface_height)).

    The function automatically names the temperatures after the depth - so you will want to rename them something sensible. You will need to pull up the original column names for your lake, and then you will need to change the column names in the rename() function accordingly!

    # pull temperature data
    LakeTemp <- get_temp(file = baseline, reference = "surface",
                                z_out = c(0, min(LakeDepth$surface_height)))
    
    # get column names
    colnames(LakeTemp)
    [1] "DateTime"              "temp_0"                "temp_33.3833669399235"
    # rename columns 
    LakeTemp <- LakeTemp %>%
      rename(Baseline_SurfaceTemp = temp_0,
             Baseline_BottomTemp = temp_33.3833669399235)
    Give it a whirl

    Visualize the change of surface and bottom temperature for the baseline simulation and briefly describe your results. Remember to be specific and to start with general patterns and then highlight important/notable details. Remember to include a title and axes labels.

    Did it!

    [Your answer here]

    Hint - to plot both temperatures in the same plot, you will need to first create a tidy data set!

    If you look at the figure below, you can see that you will need a column for date (x-axis), temperature (y-axis) and and column that contains simulated scenarios (this becomes the color coding for your line graphs).

    Once you do that, it should be straightforward to create a plot that look something like this:

    Figure 31.2: Change in surface (blue) and bottom (green) temperature [C] for baseline scenario.

    Let’s plot a heatmap of our lake thermal structure for the baseline scenario.

    plot_temp(file = baseline, fig_path = FALSE) 

    Consider this

    First, briefly explain how to interpret this heatmap (what data is code here how for the visualization?).

    Then, describe your results - include how the thermal structure changes over the course of the year. Be specific, use key terms such as thermocline, stratified, isotherm, epilomnion (above the thermocline), and hypolimnion (below the thermocline).

    Did it!

    [Your answer here]

    Next, we will extract the ice cover on the lake for our baseline simulation.

    LakeIce <- get_var(baseline, "hice") %>%
      rename(ice_baseline = hice)

    31.4 Give it a whirl

    Visualize the change of surface ice for your lake over the course of the year and then describe your results. Remember to be specific - when is there ice? How thick is it? Include a title & label your axes!

    :::

    Did it!

    [Your answer here]

    Your plot should look something like this:

    Figure 31.3: Change in ice cover [m] for baseline simulation.

    Now that we have a pretty good idea what our baseline data set looks like in terms of the thermal structure in a normal year, our next step will be to simulate what the thermal structure will look like in a typical or extreme El Nino year compared to an extreme El Nino year.

    31.5 Acknowledgments

    These activities are based on the EDDIE Teleconnections module.3

  • 3 Farrell, K.J. and C.C. Carey. 18 May 2018. Macrosystems EDDIE: Teleconnections. Macrosystems EDDIE Module 3, Version 1.